home *** CD-ROM | disk | FTP | other *** search
/ Celestin Apprentice 4 / Apprentice-Release4.iso / Source Code / Libraries / DCLAP 4j / SeqPups / apps / clustalw.src / alnscore.c next >
Encoding:
C/C++ Source or Header  |  1995-12-17  |  2.3 KB  |  118 lines  |  [TEXT/R*ch]
  1. #include <stdio.h>
  2. #include <stdlib.h>
  3. #include <math.h>
  4. #include "clustalw.h"
  5.  
  6. #define MAX(a,b) ((a)>(b)?(a):(b))
  7. #define MIN(a,b) ((a)<(b)?(a):(b))
  8.  
  9. /*
  10.  *       Prototypes
  11.  */
  12.  
  13. extern void errout();
  14. extern void *ckalloc(size_t bytes);
  15. extern void ckfree(void *);
  16. extern int  get_matrix(int *matptr, int *xref, int matrix[NUMRES][NUMRES], int neg_flag);
  17.  
  18. void aln_score(void);
  19. static int  count_gaps(int s1, int s2, int l);
  20.  
  21. /*
  22.  *       Global Variables
  23.  */
  24.  
  25. extern float gap_open;
  26. extern int   nseqs;
  27. extern int   *seqlen_array;
  28. extern int   blosum45mt[];
  29. extern int   def_aa_xref[];
  30. extern int   debug;
  31. extern int   max_aa;
  32. extern char  **seq_array;
  33.  
  34.  
  35. void aln_score(void)
  36. {
  37.   static int    *mat_xref;
  38.   static int    matrix[NUMRES][NUMRES];
  39.   static int    maxres, ngaps;
  40.   static int    l1,l2,s1,s2,c1,c2;
  41.   static int    score;
  42.  
  43.   int i;
  44.   int    *matptr;
  45.  
  46.   matptr = blosum45mt;
  47.   mat_xref = def_aa_xref;
  48.   maxres = get_matrix(matptr, mat_xref, matrix, TRUE);
  49.   if (maxres == 0)
  50.     {
  51.        fprintf(stderr,"Error: matrix blosum30 not found\n");
  52.        return;
  53.     }
  54.  
  55.   score=0.0;
  56.   for (s1=1;s1<=nseqs;s1++)
  57.    {
  58.     for (s2=1;s2<s1;s2++)
  59.       {
  60.  
  61.         l1 = seqlen_array[s1];
  62.         l2 = seqlen_array[s2];
  63.         for (i=1;i<l1 && i<l2;i++)
  64.           {
  65.             c1 = seq_array[s1][i];
  66.             c2 = seq_array[s2][i];
  67.             if ((c1>=0) && (c1<=max_aa) && (c2>=0) && (c2<=max_aa))
  68.                 score += matrix[c1][c2];
  69.           }
  70.  
  71.         ngaps = count_gaps(s1, s2, l1);
  72.  
  73.         score += 100 * gap_open * ngaps;
  74.  
  75.       }
  76.    }
  77.  
  78.   score /= 100;
  79.  
  80.   fprintf(stderr,"\nAlignment Score %d\n", (pint)score);
  81.  
  82. }
  83.  
  84. static int count_gaps(int s1, int s2, int l)
  85. {
  86.     int i, g;
  87.     int q, r, *Q, *R;
  88.  
  89.  
  90.     Q = (int *)ckalloc((l+2) * sizeof(int));
  91.     R = (int *)ckalloc((l+2) * sizeof(int));
  92.  
  93.     Q[0] = R[0] = g = 0;
  94.  
  95.     for (i=1;i<l;i++)
  96.       {
  97.          if (seq_array[s1][i] > max_aa) q = 1;
  98.          else q = 0;
  99.          if (seq_array[s2][i] > max_aa) r = 1;
  100.          else r = 0;
  101.  
  102.          if ((Q[i-1] <= R[i-1]) && (q != 0) && (1-r != 0) ||
  103.              (Q[i-1] >= R[i-1]) && (1-q != 0) && (r != 0))
  104.              g += 1;
  105.          if (q != 0) Q[i] = Q[i-1]+1;
  106.          else Q[i] = 0;
  107.  
  108.          if (r != 0) R[i] = R[i-1]+1;
  109.          else R[i] = 0;
  110.      }
  111.  
  112.    ckfree((void *)Q);
  113.    ckfree((void *)R);
  114.  
  115.    return(g);
  116. }
  117.           
  118.